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ABSTRACT 

The transitional flow over a multi-element airfoil 
in a landing configuration are computed using a two- 
equation transition model. The transition model is pre- 
dictive in the sense that the transition onset is a result 
of the calculation and no prior knowledge of the transi- 
tion location is required. The computations were per- 
formed using the INS2D Navier-Stokes code. Overset 
grids are used for the three-element airfoil. The airfoil 
operating conditions are varied for a range of angle of 
attack and for two different Reynolds numbers of 5 mil- 
lion and 9 million. The computed results are compared 
with experimental data for the surface pressure, skin 
friction, transition onset location, and velocity magni- 
tude. In general, the comparison shows a good agree- 
ment with the experimental data. 

INTRODUCTION 

Flows past multi-element airfoils have been sub- 
jected to intensive experimental and computational 
studies 1,2,3 for the past two decades. An accurate pre- 
diction of such flows can enhance the performance and 
the safety factor of aircrafts in high-lift operations. 
Even with the advances of computational fluid dynam- 
ics (CFD), the prediction of high-lift flow fields remains 
a challenge. The flow fields around multi-element air- 
foils are complex and are known to be dominated by 
different flow mechanisms at different operating condi- 
tions. 

A significant portion of the flow around multi- 
element airfoils is transitional over a wide range of op- 
erating condition. The onset location of the transition 
process varies with, for example, angles of attack and 
Reynolds numbers. In addition to boundary layer 
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transition, free shear flow transition is also likely to 
occur in, for example, the flow over the slat. As the 
slat wake convects downstream, it may further interact 
with the transitional flow over the succeeding elements, 
such as the flow over the main element and the flap. As 
a result, an accurate prediction of the flow transition 
represents a major challenge for the computation of the 
flow over multi-element airfoils and the use of a transi- 
tion model often becomes a discriminating factor in the 
outcomes of a prediction. A comprehensive evaluation 
of existing transitional empirical correlations 4 showed 
that none of the correlations evaluated gave satisfactory 
results for the rather complex high lift flows. 

Turbulence models are frequently used in the high- 
lift flow calculations. In this practice, turbulence mod- 
els are “turned-on” or “ramped-up” at a designated 
location on the airfoils. The designated locations are 
usually correlated with the measured location of tran- 
sition onset and the computational results can vary sig- 
nificantly with the assigned transition location. 

In this paper, the transitional flow fields over a 
multi-element airfoil are calculated using a k — e two- 
equation transition model 5 . The model uses an effec- 
tive eddy-viscosity by coupling an intermittence-like 
corrections to a turbulence eddy -viscosity. The inter- 
mittence function varies with the evolution of the dis- 
turbance kinetic energy, k , which is a solution variable 
for the present k — e modeling methodology. As a re- 
sult, there is no need for any prior knowledge of the 
targeted transition process. In this sense, the present 
model is a truly predictive transition model. The model 
has been applied successfully to by-pass transitional 
flows over flat plates 5 . The predicted transition onset 
and the length of the transition region were in a very 
good agreement with the measurement for all the cases 
studied. The transition model and the results of its ap- 
plication to a multi-element airfoil are presented in the 
following sections. 

The Navier-Stokes solver used in the study is de- 
scribed in the next section. The Chimera composite 
grids used are also presented. 
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FLOW SOLVER 


The flow calculations presented here were performed 
with the INS2D code 6 - 7 ’ 8 . The code solves the two 
dimensional, incompressible form of the Reynolds- 
averaged Navier-Stokes equations using an artificial 
compressibility approach. The INS2D code is a fi- 
nite difference solver, using Roe’s third order upwind- 
biased, flux-difference splitting for the convective terms 
and a second-order central differencing for the viscous 
terms. The code is capable of solving both steady and 
unsteady problems with point-to-point matched grids 
or Chimera overset grids. The details of the INS2D 
code can be found in the references cited above. The 
INS2D code has been applied extensively to the flow 
over various single- and multi-element airfoils 6,7,8 and 
has been shown to be a very reliable tool for such 
purposes. Generating the necessary grids for multi- 
element airfoil flow calculations can be very time- 
consuming regardless of the grid topology one chooses 
to use. In this study, overset grids were used and the 
grids were generated using the OVERMAGG 9 software. 
OVERMAGG is an automated script system specifi- 
cally designed to perform overset grid generation for 
multi-element airfoils. The use of OVERMAGG has 
resulted in a significant amount of saving in time dur- 
ing this study. 

The present transition model was implemented as a 
separate module to the INS2D code. There are three 
other model options in the INS2D code. These models 
have been applied to the flow over the 30P30N airfoil 
geometry which is also used in this study and the results 
have been reported 3 . 

TRANSITION MODEL 


The transitional model is developed based on the 
experimental observation that the flows in transition 
are highly intermittent. Locally, the flow is turbulent 
as the turbulent spots convect through. In between 
the passing of the turbulent spots, the flow relaxes to 
a disturbed laminar state. In this case, a measure of 
the fraction of time that the boundary layer is turbu- 
lent is the intermittency factor. Normally, the inter- 
mittent nature of the transitional flow is not accounted 
for in the development of turbulence models. Turbu- 
lence models have been in general designed for flows in 
a fully turbulent state. To be able to use a turbulence 
model in the intermittent transitional flow region, we 
incorporate a weighting factor, 7, in the turbulent eddy 
viscosity , /i*, obtained from a “parent” two-equation 


turbulence model. That is, 

thr = 1 tH 0 ) 


where fi tr denotes transitional flow eddy-viscosity. The 
7 function, or intermittency correction function, resem- 
bles the physical flow intermittency factor. Similar to 
the intermittency factor that determines the fraction 
of time the flow is turbulent in the transition region, 7 
determines the fraction of turbulent eddy-viscosity to 
be used in the transition region. Therefore, 7 should 
also vary monotonically through the transition region. 

In the present study, the intermittency correction 
function 7 is defined in terms of the streamwise varia- 
tion of the peak disturbance energy. The disturbance 
kinetic energy increases as the flow evolves from lam- 
inar to transitional state. Its local peak level also in- 
creases. 

Therefore, the disturbance energy level is an appro- 
priate parameter to characterize the flow transitional 
region. The intermittency correction function involves 
the peak disturbance energy level and the local free 
stream turbulence level. 


where 
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k e and u T denote the local disturbance kinetic energy 
in the free stream and the local frictional velocity. The 
free stream turbulence intensity was quite low in the 
wind tunnel where the transition measurement were 
made. In the present application, k e is set according to 
the turbulence level in the outer stream. The model co- 
efficient kT relates to the threshold for the formation of 
turbulent spots near the wail. It has been shown 10 that 
turbulent energetic eddies reduce to Kolmogorov eddies 
near the wall, where k+ — 0.25, and all the wall param- 
eter are characterized by Kolmogorov microscales. An 
estimate for the spot generation threshold level in terms 
of k+ is thus set at 0.25. It should be noted that the 
wall boundary condition for k was also determined by 
the Kolmogorov behavior of near-wall turbulence in the 
parent two— equation model used in this study. There- 
fore, the selection of the value of k+ is consistent with 
the current parent model, k 2 is determined by exam- 
ining fully developed turbulent boundary layer and its 
value is set at 4.5. The term k+ — k+ in Eq.(3) is used 
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to normalize the term k+ - kf , such that the transition 
progress variable, /, satisfies the equation, 


0 < / < 1. 


The effect of the local free stream turbulence is ac- 
counted for in Eqs.(4) and (5). In Eq.(4), k+ repre- 
sents the difference between the peak value and the free 
stream value of the profile at a streamwise location, 
x. In this study, the k+ profile is taken along the di- 
rection normal to the surface. Eq.(4) allows the flow 
transition to proceed according to the development of 
the internal peak level of k + . 

As was mentioned earlier, the present transition 
model coupled the intermittency formulation with a 
parent two-equation turbulence model. The parent 
model used here is a variation of the Shih and Lum- 
ley model 10 and has been used in Liou and Shih 11 and 
Liou et. al 12 . 

The model equations for k and e are, 


k t t + U{k { — [(v •+* t/ t)k > i],< — xij Uij € (6) 
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where 


Ci = 1.44, C 2 = 1.92, <r c = 1.3 
/ 2 = l-0.22exp[-(^-) 2 ], Rt= k ~ 


(7) 


C u = 


“ Ao + 
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where 

c ( *> = JSijSij + ciijCiij 

Qij is the mean rotation rate viewed in a rotating ref- 
erence frame with the angular velocity w*. The param- 
eter A s is determined by 


A, = v^6cos <j > , 

W — ^ijSjkSk, 
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The damping function is defined by 

= [1. - exp(— (difi/t 4- a 3 R 3 k + a 5 fl*))]^ 


(9) 


( 10 ) 


where 

a, = 1.5 x 10~ 3 , a 3 =10“ 9 , a 5 = 5xlO' 10 



The near-wall boundary conditions for the turbulent 
quantities are 

t jA 

jfc = 0.25u? and ^ = 0.251—^ (11) 

T v 

where u T denotes the surface friction velocity. 

Note that no measurement was made to determine 
the turbulent intensity and the decay of turbulent ki- 
netic energy in the experiment where the present data 
were taken. These information are needed to determine 
the value of the disturbance kinetic energy and its dis- 
sipation rate at the inlet of the computational domain 
in the present approach. A turbulent intensity level of 
0.05% was cited and used in the present analysis and 
the equivalent eddy— viscosity was set at the laminar 
level. 

The present intermittency correction function tech- 
nique uses the variation of the disturbance energy to 
computationally characterize the transition progress. 
Therefore, the calculated intermittency correction func- 
tion does not necessarily correspond to the physical flow 
intermittency factor. 

GEOMETRY 


The McDonnell Douglas 30P-30N landing config- 
uration was used in this study. Figure 1 shows the 
geometry and the stations where profile data will be 
presented. Both the leading-edge slat and the trailing- 
edge flap have a deflection angle of 30°. The airfoil 
has been tested extensively at NASA Langley Low- 
Turbulence Pressure Tunnel (LTPT) 1,13 . The geometry 
has been used as a test case in a CFD Challenge Work- 
shop held at NASA Langley in 1993 and many compu- 
tational results obtained by using different solvers and 
models have been reported 2,3,8 . 

The computational grid consists of seven zones and 
a total of about 103,000 grid points. Figure 2 shows 
the C grids around the slat, the main element, and the 
flap, which consist of 175x45, 38 lx 133, and 221x85 
points, respectively. Two H grids were used for the 
main-element cove and for the wake region of the flap. 
A background H grid extends to the top and the bot- 
tom walls of the wind tunnel. When the grid for one 
element overlaps another element (s), holes and outer 
boundaries are determined by using the PEGSUS code. 
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The communication between the overset grids are also 
handled by the PEGSUS code. 

Comprehensive grid independence studies have 
been performed and reported previously using similar 
overset grids for the same geometry 3 8 These and other 
calculations for the same geometry indicate that the 
grid point distributions in the wake regions can have 
certain effects on the velocity profiles in those regions 
and the wall pressure is not sensitive to the grid den- 
sity. These results also show that further refinement 
of the grid from the current level does not change the 
results appreciably. 

In the following, results will be presented for two 
Reynolds numbers of five million and nine million and 
a range of angle of attack (AOA) for each Reynolds 
number. 

RESULTS AND DISCUSSIONS 

Figure 3 shows a comparison of the surface pressure 
coefficient, C p , distributions over the slat, main ele- 
ment, and flap for Re=9 million and AOA=19°. The 
calculated distributions agree well with the measure- 
ment. The surface pressure coefficients for A0A=8° 
are presented in Figure 4. The predicted pressure over 
the suction side of the the slat is slightly lower than the 
data. The overall agreement of the predicted surface 
pressure with the measurement is satisfactory for both 
cases, showing that the transitional model is capturing 
the outer “inviscid” flow fields. 

Figure 5 shows the calculated and the measured skin 
friction coefficients over the main element for Re=9 mil- 
lion and AOA=19°. The computed results agree very 
well with the available measured data. 

The flow around the slat is complicated and is diffi- 
cult to measure due to the physical size of the slat. The 
transition of the boundary layer from laminar to tur- 
bulent flow on the upper suction surface has long been 
recognized as being one of the critical flow features that 
are yet to be modeled correctly. In the experiment the 
onset and the extend of a transition process were deter- 
mined based upon the skewness, flatness, and standard 
deviation of hot film data. The transition parameters 
indicated by the three different measured signals are 
not always the same and are subject to interpretation. 
Also, these variables are not derivable from the results 
of the present computational methodology. 

In this study, the onset location of flow transition 
is assumed to be indicated by the minimum value of 
the calculated Cj. Figure 6 shows the calculated C / 
over the slat for Re=9 million and AOA=19°. The Cj 
obtained by using the parent turbulence model is also 


included for comparison. 

Figure 7 shows a comparison of the predicted and 
the measured transition onset locations on the suction 
side of the three elements for Re=9 million and AOA= 
8°, 10°, 12°, 16°, and 19°. The present model gives 
rather satisfactory predictions of the transition loca- 
tions on all three elements for the five AOA measured 
at this Reynolds number. 

The variations of the transition onset locations with 
AOA on the suction surface of the individual element 
are shown in Figure 8. Figure 8. a shows the slat results. 
The range of x/c shown corresponds to the slat chord 
length. The measured data indicated a change of tran- 
sition location toward the leading edge of the slat as the 
AOA increases. It is apparent that the present model 
has predicted this trend well. Moreover, the predicted 
transition locations also agree well with the measured 
data. Figures 8.b and 8.c show comparisons between 
the predicted and the measured transition locations for 
the main element and flap, respectively. Similar to Fig- 
ure 8. a, the range of x/c shown correspond to their re- 
spective chord lengths. The predicted locations of tran- 
sition onset on the main element and flap move little 
with the change of AOA, which agree with the measure- 
ment. The predicted onsets occur slightly upstream of 
the measured locations for both the main and flap. 

Another view of the transition locations is given in 
Figure 9, which shows the locations of the predicted and 
the measured transition onset on the 30P30N airfoil for 
AOA=19°. 

Figure 10 shows a comparison of the predicted and 
the measured transition onset locations on the suction 
side of the three elements for Re=5 million and AOA= 
8°, 10°, 12°, 16°, and 19°. For this lower Reynolds 
number, the present model also predicts onset locations 
that are in a good agreement with the measured data 
for all three elements. 

Figure 11 shows the variation of the transition onset 
location with AOA for the individual element. Except 
for the lowest AOA of 8°, the present model predicts 
reasonably well the transition locations on the suction 
side of the slat. The flap transition locations were also 
satisfactorily predicted. 

Note that the present transition model predicts the 
transition location, rather than using it as an input. 
Through the ranges of angle of attack and Reynolds 
number tested, the present model has predicted satis- 
factorily the trend of change of the transition location 
with angle of attack on all three elements of the 30P30N 
geometry. The absolute locations of transition onset are 
also predicted in a fairly consistent manner. 

The present transition model uses the evolution of 
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the internal peak of the calculated disturbance kinetic 
energy and the variation of the disturbance in the outer 
stream to characterize the development of a transitional 
flow. For a multi-element airfoil, such as the 30P30N 
airfoil, the wakes of the preceding elements convect 
downstream and further interact with the transitional 
flow over the succeeding elements. Figure 12 shows 
the calculated distributions of the disturbance kinetic 
energy, k , along the direction normal to the upper sur- 
face of the slat, main, and flap element, respectively, 
for Re=9 million and A0A=19°. The local high val- 
ues in the outer stream represent the high disturbance 
energy associated with the slat and main wakes. The 
wakes generated by the preceding elements can also be 
clearly observed in the velocity magnitude contour plot 
shown in Figure 13. For this operating condition, the 
slat wake can be identified as far downstream as in the 
flap region with little diffusion. It can be argued that 
there are receptive mechanisms for the boundary layer 
on the main, or flap, element to be perturbed by the 
high-intensity disturbances in the slat, or main, wake. 
In this study, the effect of the wake generated by the 
preceding element(s) is taken into account by using the 
highest k in the wake(s) as the k e in Eq.(4) for the suc- 
ceeding element(s). Near the leading edge of the main 
element, where the surface normal intersects with the 
slat cove, k e can vary significantly due to the separated 
shear layer in the cove region. While not affecting the 
wall pressure, it causes low amplitude variation in the 
skin friction and greatly delay convergence locally. In 
this application, an average value of k e over the inter- 
sect region was used for both the main element and the 
flap. 

Figure 14 shows the velocity profiles on the main el- 
ement and on the flap for Rc=9 million and A0A=19°. 
At x/c = 0.1075, there is almost an uniform off- 
set between the data and the predictions. This 
may have been caused by possible improper data 
calibration 2 . The predicted slat wake is located near 
the measured wake immediately behind the slat trail- 
ing edge (x/c=0. 1075). The calculated values for the 
wake deficit agree reasonably well with the data at 
z/c=0.1075 and 0.45. The measured wake is wider at 
x/c=0.45. The merging of the calculated slat wake and 
the main boundary layer is apparently yet to occur at 
x/c = 0.85. The boundary layer velocity profiles are 
predicted well by the present transition model, which 
results in a better prediction of the main wake imme- 
diately behind the main element trailing edge (Figure 
14. d, x/c — 0.89817). 

Over the flap (Figures 14. e and 14. f, x/c — 1.0321 
and 1.1125), the predicted development of the slat wake 


is satisfactory. The predicted deficit of the main ele- 
ment wake is higher than the measured values. The 
width of the calculated main element wake agrees rea- 
sonably well with the data. 

Figure 15 shows the velocity profiles on the main 
element for Re=9 million and AOA=8°. The predicted 
main element surface boundary layer is more developed 
than what the data has shown at ar/c=0. 10765. There 
is a good agreement between the measurement and the 
calculated velocity profiles at x/c=0.45 and 0.85, where 
both the prediction and the measurement show a con- 
fluent wake/boundary layer. 

Overall, the predicted velocity profiles agree reason- 
ably well with the measurement. 

CONCLUDING REMARKS 

The transitional flows over a two-dimensional 
multi-element airfoil have been computed using a pre- 
dictive transition model. The model couples an inter- 
mittency function correction to a parent two-equation 
turbulence model through an effective eddy-viscosity 
for the transitional flow. Comparisons with available 
measured data for two Reynolds numbers and a range 
of angle of attack were made. These include the sur- 
face pressure, skin friction, and velocity profile. The 
predicted transition onset on all the elements were also 
compared with the experimental data. The compar- 
isons show that the present transition model predicts 
fairly well the locations of transition onset for the cases 
tested. 
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